source("utils.R")
library(DT)
library(reshape2)
library(ggplot2)
library(vegan)
library(gganimate)
library(gifski)
library(gridExtra)| AB | Class |
|---|---|
| Ciprofloxacin | Fluoroquinolone |
| Tetracyclin | Polyketide |
| Vancomycin | Glycopeptide |
| ID | phylum | species |
|---|---|---|
| YL44 | Verrucomicrobia | A. muciniphila |
| I48 | Bacteroidetes | B. caecimuris |
| YL27 | Bacteroidetes | M. intestinale |
| YL45 | Proteobacteria | T. muris |
| YL2 | Actinobacteria | B. longum |
| KB1 | Firmicutes | E. faecalis |
| KB18 | Firmicutes | A. muris |
| YL32 | Firmicutes | C. clostridioforme |
| YL31 | Firmicutes | F. plautii |
| YL58 | Firmicutes | B. coccoides |
| I49 | Firmicutes | L. reuteri |
| I46 | Firmicutes | C. innocuum |
Define colors for species
species_cols <- c("KB1" = "#73F440",
"YL2" = "#1F3C8E",
"KB18" = "#616160",
"YL27" = "#F49D44",
"YL31" = "#42933A",
"YL32" = "#3D98B1",
"YL44" = "#AC80B6",
"YL45" = "#E43531",
"I46" = "#000100",
"I48" = "#C0622C",
"I49" = "#4FAF75",
"YL58" = "#193E1E")
species_cols_names <- c("E.faecalis" = "#73F440",
"B.longum" = "#1F3C8E",
"A.muris" = "#616160",
"M.intestinale" = "#F49D44",
"F.plautii" = "#42933A",
"C.clostridioforme" = "#3D98B1",
"A.muciniphila" = "#AC80B6",
"T.muris" = "#E43531",
"C.innocuum" = "#000100",
"B.caecimuris" = "#C0622C",
"L.reuteri" = "#4FAF75",
"B.coccoides" = "#193E1E")
species_names <- read.table("members.txt", sep = ";", header = T)Raw data are copies/g feces is located in values.csv. Detection limits are included and located in lod.csv.
# read qPCR data
dat <- read.table("values.csv", header = T, sep = ";")
dat$universal <- NULL
dat.m <- melt(dat, id.vars = c("day", "mouse2" ,"group", "mouse", "applicaiton"))
# replace species ID with full nmae
dat.m$variable <- species_names[match(dat.m$variable, species_names$ID),]$species
# load detection limit
lod <- read.table("lod.csv", sep = ";", header = T)
lod <- melt(lod, id.vars = "type")
# change id to species name
lod$variable <- species_names[match(lod$variable, species_names$ID),]$species
lod <- lod[which(lod$variable != "universal"),]
lod_upper <- lod[which(lod$type == "upper"),]
lod_lower <- lod[which(lod$type == "lower"),]
lod2 <- read.table("lod.csv", header = T, sep = ";")
rownames(lod2) <- c("upper", "lower")
lod2$type <- NULL
lod2 <- as.data.frame(t(lod2))dat.summary <- aggregate(. ~ day, mean, data = dat.m)
treatment <- data.frame(begin = c(1, 15),
end = c(10, 20))
p <- ggplot(dat.m, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_grid(variable ~ group) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
pdat.m_subset <- dat.m[which(dat.m$group == "water"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_wrap(~variable) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("absolute abundance in control group")
pdat.m_subset <- dat.m[which(dat.m$group == "Ciprofloxacin"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_wrap(~variable) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("absolute abundance in Ciprofloxacin group")
pdat.m_subset <- dat.m[which(dat.m$group == "Tetracyclin"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_wrap(~variable) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("absolute abundance in Tetracyclin group")
pdat.m_subset <- dat.m[which(dat.m$group == "Tetracyclin"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
#p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
#p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_wrap(~variable, scales = "free") + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("absolute abundance in Tetracyclin group")
pdat.m_subset <- dat.m[which(dat.m$group == "Vancomycin"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + facet_wrap(~variable) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("absolute abundance in Vancomycin group")
pdat.m_subset <- dat.m[which(dat.m$group == "Tetracyclin"),]
dat.m_subset <- dat.m_subset[which(dat.m_subset$variable == "M.intestinale"),]
p <- ggplot(dat.m_subset, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
#p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
#p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p <- p + ggtitle("C.innocuum in Tetracyclin group")
p <- p + ylim(c(min(dat.m_subset$value), max(dat.m_subset$value)))
pdat <- read.table("values.csv", header = T, sep = ";")
dat$universal <- NULL
rownames(dat) <- paste0(dat$group, "-", dat$mouse2, "-", dat$day)
type <- dat$group
day <- dat$day
dat[,c(1:5)] <- NULL
day2 <- rep("none", length(day))
day2[which(day == 0)] <- "before treatment"
day2[which(day == 4)] <- "treatment response"
day2[which(day == 18)] <- "treatment response"
day2[which(day == 53)] <- "treatment response"
day2[which(day == 67)] <- "treatment response"
day2[which(day > 67)] <- "after treatment"
day2[which(day2 == "none")] <- "resilience"
dat_t <- t(as.data.frame(dat))
dat_scaled <- scale(dat_t, center = FALSE, scale = colSums(dat_t))
dat <- t(as.data.frame(dat_scaled))
# generating distance matrix
dist <- vegan::vegdist(dat, method = 'bray')
pcoa <- cmdscale(dist, k = 2, eig = TRUE)
# generate the data frame for plotting, consisting of coordinates and the annotation
coordinates <- data.frame(
pcoa1 = pcoa$points[,1],
pcoa2 = pcoa$points[,2],
type = type,
round = day2,
time = day,
row.names = NULL
)
# calculated explained variance
eigenvalues <- eigenvals(pcoa)
variance <- eigenvalues / sum(eigenvalues)
variance1 <- 100 * round(variance[1], digits = 3)
variance2 <- 100 * round(variance[2], digits = 3)
# plotting
p <- ggplot(coordinates, aes(x = pcoa1, y = pcoa2, color = day2))
p <- p + geom_point(alpha = 1, stroke = 0, aes(shape = type), size = 1)
#p <- p + scale_size(range = c(2, 6), name = "time (day)")
p <- p + xlab(paste("PCo 1 (", round(variance1), "%)", sep = ""))
p <- p + ylab(paste("PCo 2 (", round(variance2), "%)", sep = ""))
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + scale_colour_brewer(palette = "Set1", name = "state", direction = -1)
p <- p + scale_shape_manual(values = c(16, 17, 15, 1), name = "group")
p <- p + theme_bw() + geom_hline(yintercept = 0, linetype = 3)
p <- p + geom_vline(xintercept = 0, linetype = 3)
p <- p + theme(aspect.ratio = 1)
pdat <- read.table("values.csv", header = T, sep = ";")
dat$universal <- NULL
rownames(dat) <- paste0(dat$group, "-", dat$mouse2, "-", dat$day)
type <- dat$group
day <- dat$day
dat[,c(1:5)] <- NULL
day2 <- rep("none", length(day))
day2[which(day == 0)] <- "before treatment"
day2[which(day == 4)] <- "treatment response"
day2[which(day == 18)] <- "treatment response"
day2[which(day == 53)] <- "treatment response"
day2[which(day == 67)] <- "treatment response"
day2[which(day > 67)] <- "after treatment"
day2[which(day2 == "none")] <- "resilience"
# generating distance matrix
dist <- vegan::vegdist(dat, method = 'bray')
pcoa <- cmdscale(dist, k = 2, eig = TRUE)
# generate the data frame for plotting, consisting of coordinates and the annotation
coordinates <- data.frame(
pcoa1 = pcoa$points[,1],
pcoa2 = pcoa$points[,2],
type = type,
round = day2,
time = day,
row.names = NULL
)
# calculated explained variance
eigenvalues <- eigenvals(pcoa)
variance <- eigenvalues / sum(eigenvalues)
variance1 <- 100 * round(variance[1], digits = 3)
variance2 <- 100 * round(variance[2], digits = 3)
# plotting
p <- ggplot(coordinates, aes(x = pcoa1, y = pcoa2, color = day2))
p <- p + geom_point(alpha = 1, stroke = 0, aes(shape = type), size = 1)
#p <- p + scale_size(range = c(2, 6), name = "time (day)")
p <- p + xlab(paste("PCo 1 (", round(variance1), "%)", sep = ""))
p <- p + ylab(paste("PCo 2 (", round(variance2), "%)", sep = ""))
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + scale_colour_brewer(palette = "Set1", name = "state", direction = -1)
p <- p + scale_shape_manual(values = c(16, 17, 15, 1), name = "group")
p <- p + theme_bw() + geom_hline(yintercept = 0, linetype = 3)
p <- p + geom_vline(xintercept = 0, linetype = 3)
p <- p + theme(aspect.ratio = 1)
pggplot(coordinates, aes(x = pcoa1, y = pcoa2, colour = type, size = day, shape = type)) +
geom_point(alpha = 0.7, show.legend = TRUE) +
scale_colour_brewer(palette = "Set1", name = "group", direction = 1) +
theme_bw() + geom_hline(yintercept = 0, linetype = 3) +
scale_shape_manual(values = c(16, 17, 15, 1), name = "group") +
geom_vline(xintercept = 0, linetype = 3) + scale_size(range = c(2, 6)) +
theme(aspect.ratio = 1) +
# Here comes the gganimate specific bits
labs(title = 'time (days): {frame_time}', x = 'PCo 1', y = 'PCo 2') +
transition_time(time) +
ease_aes('linear') +
enter_fade() +
exit_fade()dat <- read.table("values.csv", header = T, sep = ";")
dat$universal <- NULL
rownames(dat) <- paste0(dat$group, "-", dat$mouse2, "-", dat$day)
type <- dat$group
day <- dat$day
dat[,c(1:5)] <- NULL
day2 <- rep("none", length(day))
day2[which(day == 0)] <- "before treatment"
day2[which(day == 4)] <- "treatment response"
day2[which(day == 18)] <- "treatment response"
day2[which(day == 53)] <- "treatment response"
day2[which(day == 67)] <- "treatment response"
day2[which(day > 67)] <- "after treatment"
day2[which(day2 == "none")] <- "resilience"
dat_t <- t(as.data.frame(dat))
dat_scaled <- scale(dat_t, center = FALSE, scale = colSums(dat_t))
dat <- t(as.data.frame(dat_scaled))
# generating distance matrix
dist <- vegan::vegdist(dat, method = 'bray')
pcoa <- cmdscale(dist, k = 2, eig = TRUE)
# generate the data frame for plotting, consisting of coordinates and the annotation
coordinates <- data.frame(
pcoa1 = pcoa$points[,1],
pcoa2 = pcoa$points[,2],
type = type,
round = day2,
time = day,
row.names = NULL
)
ggplot(coordinates, aes(x = pcoa1, y = pcoa2, colour = type, size = day, shape = type)) +
geom_point(alpha = 0.7, show.legend = TRUE) +
scale_colour_brewer(palette = "Set1", name = "group", direction = 1) +
theme_bw() + geom_hline(yintercept = 0, linetype = 3) +
scale_shape_manual(values = c(16, 17, 15, 1), name = "group") +
geom_vline(xintercept = 0, linetype = 3) + scale_size(range = c(2, 6)) +
theme(aspect.ratio = 1) +
# Here comes the gganimate specific bits
labs(title = 'time (days): {frame_time}', x = 'PCo 1', y = 'PCo 2') +
transition_time(time) +
ease_aes('linear') +
enter_fade() +
exit_fade()dat <- read.table("values.csv", header = T, sep = ";")
dat$universal <- NULL
rownames(dat) <- paste0(dat$group, "-", dat$mouse2, "-", dat$day)
type <- dat$group
day <- dat$day
dat[,c(1:5)] <- NULL
day2 <- rep("none", length(day))
day2[which(day == 0)] <- "before treatment"
day2[which(day == 4)] <- "treatment response"
day2[which(day == 18)] <- "treatment response"
day2[which(day == 53)] <- "treatment response"
day2[which(day == 67)] <- "treatment response"
day2[which(day > 67)] <- "after treatment"
day2[which(day2 == "none")] <- "resilience"
dat_t <- t(as.data.frame(dat))
dat_scaled <- scale(dat_t, center = FALSE, scale = colSums(dat_t))
dat <- t(as.data.frame(dat_scaled))
# generating distance matrix
dist <- vegan::vegdist(dat, method = 'bray')
pcoa <- cmdscale(dist, k = 2, eig = TRUE)
# generate the data frame for plotting, consisting of coordinates and the annotation
coordinates <- data.frame(
pcoa1 = pcoa$points[,1],
pcoa2 = pcoa$points[,2],
type = type,
round = day2,
time = day,
row.names = NULL
)
generateOnePCOA <- function(dat, day=4) {
dat_day <- coordinates[which(coordinates$time == as.numeric(day)),]
p <- ggplot(dat_day, aes(x = pcoa1, y = pcoa2, colour = type, shape = type))
p <- p + geom_point(alpha = 0.7, show.legend = TRUE)
p <- p + xlim(min(dat$pcoa1), max(dat$pcoa1))
p <- p + ylim(min(dat$pcoa2), max(dat$pcoa2))
p <- p + scale_colour_brewer(palette = "Set1", name = "group", direction = 1)
p <- p + theme_pmuench() + geom_hline(yintercept = 0, linetype = 3)
p <- p + scale_shape_manual(values = c(16, 17, 15, 1), name = "group")
p <- p + geom_vline(xintercept = 0, linetype = 3) + scale_size(range = c(2, 6))
p <- p + theme(aspect.ratio = 1) + labs(title = paste0('time (days): ', as.character(day)), x = 'PCo 1', y = 'PCo 2') + theme(legend.position = "none")
return(p)
}
d0 <- generateOnePCOA(coordinates, day = 0)
d4 <- generateOnePCOA(coordinates, day = 4)
d44 <- generateOnePCOA(coordinates, day = 44)
d67 <- generateOnePCOA(coordinates, day = 67)
pdf(file="pcoa_time.pdf")
grid.arrange(d0, d4, d44, d67, nrow = 2)
dev.off()## quartz_off_screen
## 2
div <- as.data.frame(diversity(dat))
colnames(div) <- "diversity"
div$group <- type
div$round <- day2
data_summary <- function(x) {
m <- mean(x)
ymin <- m - sd(x)
ymax <- m + sd(x)
return(c(y = m,ymin = ymin,ymax = ymax))
}
# draw
p <- ggplot(div, aes(x = day, y = diversity, color = group))
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_jitter(size = 1) # + geom_violin(trim = FALSE)
p <- p + scale_shape_manual(values = c(16, 17, 15, 1), name = "group")
p <- p + theme_pmuench() + xlab("time (days)") + ylab("shannon diversity")
p <- p + scale_colour_brewer(palette = "Set1", name = "group", direction = 1)
p <- p + stat_summary(fun.data = data_summary, color = "black")
pdat <- read.table("values.csv", header = T, sep = ";")
mouse <- dat$mouse2
# normalize by 16 rRNA gene copy number
dat$I46 <- dat$I46 / 5
dat$I48 <- dat$I48 / 3
dat$I49 <- dat$I49 / 6
dat$KB1 <- dat$KB1 / 4
dat$KB18 <- dat$KB18 / 2
dat$YL2 <- dat$YL2 / 3
dat$YL27 <- dat$YL27 / 4
dat$YL31 <- dat$YL31 / 2
dat$YL32 <- dat$YL32 / 5
dat$YL44 <- dat$YL44 / 3
dat$YL45 <- dat$YL45 / 6
dat$YL58 <- dat$YL58 / 5
dat$universal <- NULL
dat[,c(1:5)] <- NULL
dat_t <- t(as.data.frame(dat))
dat_scaled <- scale(dat_t, center = FALSE, scale = colSums(dat_t))
dat <- as.data.frame(t(as.data.frame(dat_scaled)))
dat$group <- type
dat$day <- day
dat$mouse <- mouse
df.m <- melt(dat, id.vars = c("day", "group", "mouse"))
df.m$value <- as.numeric(df.m$value)
p <- ggplot(df.m, aes(x = day, y = value, fill = variable))
p <- p + facet_grid(mouse ~ group)
#p <- p + geom_bar(stat = "identity")
p <- p + geom_area()
p <- p + geom_rect(data = df.m, xmin = 0, xmax = 4,ymin = 0, ymax = 1, color = "grey50", fill = NA, alpha = .5, size = 0.2)
p <- p + geom_rect(data = df.m, xmin = 0, xmax = 4,ymin = 0, ymax = 1, color = "grey50", fill = NA, alpha = .5, size = 0.2)
p <- p + geom_rect(data = df.m, xmin = 14, xmax = 18,ymin = 0, ymax = 1, color = "grey50", fill = NA, alpha = .5, size = 0.2)
p <- p + geom_rect(data = df.m, xmin = 49, xmax = 53,ymin = 0, ymax = 1, color = "grey50", fill = NA, alpha = .5, size = 0.2)
p <- p + geom_rect(data = df.m, xmin = 64, xmax = 68,ymin = 0, ymax = 1, color = "grey50", fill = NA, alpha = .5, size = 0.2)
p <- p + theme_pmuench()## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
# read qPCR data
dat <- read.table("values.csv", header = T, sep = ";")
dat <- dat[,-c(6:17)] # remove non-universal
dat.m <- melt(dat, id.vars = c("day", "mouse2" ,"group", "mouse", "applicaiton"))
p <- ggplot(dat.m, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
p <- p + facet_wrap(~group, nrow = 2) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p# read qPCR data
dat <- read.table("values.csv", header = T, sep = ";")
dat.m <- melt(dat, id.vars = c("day", "mouse2" ,"group", "mouse", "applicaiton"))
dat.m$variable <- species_names[match(dat.m$variable, species_names$ID),]$species
p <- ggplot(dat.m, aes(x = day, value, color = mouse2))
# annotate AB application
p <- p + geom_rect(xmin = 0, xmax = 4,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 14, xmax = 18,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 49, xmax = 53,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
p <- p + geom_rect(xmin = 64, xmax = 68,ymin = -Inf, ymax = Inf, color = "grey90", fill = "grey90", alpha = .1)
# annotate LOD
p <- p + geom_hline(data = lod_upper, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + geom_hline(data = lod_lower, aes(yintercept = value), linetype = 1, size = 0.2, color = "black")
p <- p + stat_summary(fun.y = mean, fun.ymin = mean, fun.ymax = mean, geom = "line", color = "black")
p <- p + geom_jitter(size = 0.5, show.legend = FALSE)
p <- p + scale_colour_brewer(palette = "Set1") + theme_pmuench()## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
p <- p + facet_grid(variable ~ group) + scale_y_log10()
p <- p + xlab("time (days)") + ylab("copies/g feces (log10)") #+ annotation_logticks( sides = "l")
p <- p + theme( panel.background = element_blank(),panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),plot.background = element_blank())# + geom_smooth()
p <- p + theme(legend.title = element_blank()) +
theme(strip.text.y = element_text(size = 6, colour = "black", face = "italic",angle = 0))
p